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When generalizing recent various quantum mechanical models of heavy quarkonia to Quantum 
Filed theoretical approach based on Bethe-Salpeter equation one is faced to the solutions that 
do not exist in nonrelativistic limit. Mainly, there is unexpected doubling of the spectrum when 
comparing to the experimentally known spectrum as well as the ones obtained from the solution 
of the Schroedinger equation. These additional states are not apriory unphysical as both of them 
have the same symmetry. Our study strongly suggests that these solutions appear due to the 
sensitivity of BSE to the details of the analytical form of the constituents quark and antiquark 
propagators, more specifically they are consequence of using unconfining free propagators. To show 
this explicitly we develop and describe the efficient method of the numerical solution for quarkonium 
BSE and numerically solve it for the case of pseudoscalar charmonia. For the bare propagators of 
constituents we are able to find BSE solution for arbitrarily high excited state without any 3- 
dimensional reduction. Unlike to the Schroedinger equation the excited states are not orthogonal. 
Using free charm quark propagators we observe that the ground state rjc is situated slightly above 
naive quark-antiquark threshold 2mc, thus the all excited states are situated above this threshold. 
In the second part of the paper we consider the model of BSE with confined propagators and show 
the influence on the spectrum directly in the Minkowski space for the first time. 
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1^ I. INTRODUCTION 

Cli 

(~| , Excited meson spectroscopy is a keystone experimental output, which is essential for understanding of quark- 
antiquark interaction. The dynamics of its constituents is dominantly driven by solely known strongly interacting 
quantum field theory- quantum chromodynamics. The confining forces are responsible for the absence of colored 
hadron constituents and simultaneously they lead to the spectra of the angularly and radially excited resonances. 
A lot is known from heavy quarkonia production in e^e~ anihilations, where vector quarkonia are straightforwardly 
PsJ . produced. BABAR ^T^, Belle and BESS experiments continue collection of various meson experimental data. The 

■ masses of up to day observed excited charmonia (bottomonia) lie in the large range 3 — 4.5GeF(9.5 — WGeV), which 
is a clear signal of the importance of relativistic description. One expects the similar for other spin mesons, however 

[ — ■ they are not easily accessible by the experiments. The heaviest candidate for excited pseudoscalar rjc meson was 
', X(3940) observed by Belle @ in double charmonium production is very likely f?(3S') state. 

■ Recently in CLEO Q and Belle experiments progress in the experimental determination of electromagnetic tran- 
" _ [ \ sitions has provided new data, which should be confronted with fully Poincare invariant description. Recall in naive 

^ . quantum mechanical quarkonium picture tplnS) — >■ 77(71 S)"f transitions should vanishes in zero recoil 7 for different 
" n,n due to the orthogonality of radial wave functions (see for instance [1, Q). This is not the case of relativistic 
treatment in BSE framework, where this so called "hindered" transitions are not expected to be vanishing even for 



X 



■ small photon energy. 



Historically, following ideas of Ref. confinement in heavy flavor hadron sector has typically been associated 

with a linearly or logarithmically rising potential between constituents [l^ [isj . The spin degeneracy observed in 
meson spectrum tell us that the confining part of the interaction should be largely spin independent- it should be 
Lorentz scalar. An Euclidean approximation of lattice calculations leaded to various predictions for the potential 
between heavy quark-antiquark states. While the absence of free quarks and gluons is an experimental fact and 
should be independent on the choice of gauge used in the practice, this is not the case of the potential. For instance 
Coulomb gauge potential rises linearly with the slope larger then required by Regge trajectory of light mesons, whilst 
all interpolating gauges between Landau and Coulomb gauge produce potential which is bounded from above p3 |. 
Notably, it is most flat (and naively less confining) in the Landau gauge. Within a lattice predicted potentials the 
meson spectra were commonly calculated by using relativistic or nonrelativistic quantum mechanics. Recall, relatively 
small hyperfine and fine-structure splitting in quarkonium levels is due to the Lorentz vector part of the interquark 



interaction [15l - ll8{ . which should be added for more accurate description. 
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In quantum mechanic all the energy of excited meson must be situated bellow the maximum of quantum mechanical 
potential and the states belonging to various meson masses are mutually orthogonal. On the other side, using quantum 
mechanical picture itself one can estimate validity of nonrelativistic description. Within the model of logarithmic 
confinement one has < >= 0.24 for quark of mass rric = 1.5GeV bounded inside of IS' charmonia. 

Thus a nonrelativistic description for charmonium is quite crude (this is less urgent for bottomonium where one has 
< >= 0.08). 

In Quantum Field Theory the two body bound state is described by the three-point bound state vertex function 
or, equivalently, by the Bethe-Salpeter (BS) amplitudes. Both of them are solutions of the corresponding covariant 
four-dimensional Bethe-Salpeter equation (BSE). In this framework the excited states of meson are not described by 
orthogonal wave functions, albeit the structure of BSA has definite structure dictated by the total spin of the meson. 
Furthermore, one can expect the states with properties that do not appear in the nonrelativistic limit (recall here, 
in principle th nonrelativistic reduction is not defined without certain ambiguities, since in Quantum Field Theory 
two body state is necessarily described by the equation including two time variables instead of one universal time 
in the Schroedinger equation) . The main motivations of presented study is not to provide the best fitting model of 
the mesons, but to fill missing gap in the literature and actually find the model based on the full 4-dimensional BSE 
which as far as possible realistically describes the quarkonia observed in the experiments. While not easily visible 
when study the ground states ground states, when approaching excited states, then one is faced to a large sensitivity 
of BSE with respect to the analytical structure of its kernel and the propagator of constituents as well. In this paper 
we will show one of such peculiar phenomena- the population doubling of radially excited states. 

On the other side , one can largely exploit a cumulated knowledge as many various three-dimensional reductions of 
BSE or Schroedinger equation have been often employed to extract the spectra and wave functions. Recently in Refs. 
[23 - l27| it has been found that the meson spectroscopy is better described by " confining" potential which is bounded 
from above. Irrespective of the gauge, the flatness of linearly rising potential arises form the string breaking scenario: 
including light quarks into the game then the creation of the quark- antiquark pairs, i.e. pions and other light mesons 
is energetically favorable. Such screening is in agreement with observed high radially excited heavy mesons- without 
any doubt, the spectrum deviates from the linear Regge trajectory. 

In the recent paper [28i] vector charmonia has been considered in the approximation with one dominant component. 
In presented paper we continue study for pseudoscalar cc mesons by using the BSE model. We confirm that the 
known dominant component for the ground state reminds leading for all studied excited states. To this point, 
the BSE for all components completely included has been solved for arbitrary excited state for the first time. Of 
course, comparing to vectors, the experimental knowledge of pseudoscalar charmonia is quite pure due to the worse 
experimental accessibility, thus the states with masses above 3.9GeV are a prediction of our presented BSE model. 

The second part of the presented paper is devoted to the novel treatment of BSE for mesons. The main purpose of 
this section is to show the origin of above mentioned doubling phenomena, however the model is interesting in its own. 
We use phenomenological form of confined quark propagators, which prevents the colored quarks be on mass-shell. We 
assume the quarks are short living colored excitations characterized by the propagators with two complex conjugated 
poles (or the sum of). For this purpose we accommodate a special form of the propagator proposed in the papers 
[42,,, ,43. ] , where the idea was developed for gluons. The model presented here has an amazing property: due to the 
absence of the real poles in the BSE kernel it becomes numerically solvable directly in the Minkowski space. For 
the first time in practice, when solving BSE for mesons with confining propagators, we actually avoid the Euclidean 
metric by staying in the original Minkowski space. 

II. BSE FOR PSEUDOSCALAR CHARMONIUM 

The model of BSE here has been in fact already considered for the calculation of vector charmonia, however here 
we are dealing with pseudoscalars instead. In this paper we prefer to work with BSE vertex function rather then BSE 
wave function. The BSE reads 

r(g, ^) = / iSiq P/mP, 'l)S{q + P/2)r Vik, q, P^k.i,, , (1) 

where Latin letters represent Dirac indices. Explicitly for pseudoscalar we have 

rp(g, P) = 75 {A{q, P)+ ^C{q, P)+ /iBiq, P) + [^, ^]D{q, P)) , (2) 
Alternatively the BSE vertex function can be replaced by (more singular) BS wave function x through its definition 

V{q,P) = S~\q-P/2)x{q,P)S-\q + P/2) , (3) 
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which is usually considered when performing nonrelativistic limit. The interaction kernel is given by the infinite sum 
of two-particle irreducible quark-antiquark scattering graphs in color-anticolor channel. According to the previous 
discussion the dominant part should be a Lorentz scalar, further vectorial interaction is naturally assumed in QCD. 
The interaction phenomenologically chosen here thus reads 



v., = 



(g2 - ;,2)2 



(4) 



This scalar part of the kernel is Lorentz invariant generalization of the "screened linear" potential, since it gives 
negative exponential potential in the nonrelativistic limit, Vs{r) ~ —e^^^". Such screened potentials have been 
proposed and successfully describe heavy quarkonia and heavy flavor mesons, as well as light hadron (2^ . 

The color Coulomb potential is usually taken due to the exchange of gluon. According to findings and recent 
knowledge in nonperturbative QCD the vectorial interaction should reflect two findings: the free zing of effective 
coupling in the infrared and gluon propagator suppression due to the soft gluon mass generation |3ll433j through the 
Yang-Mills Schwinger mechanism. Here we simplify and approximate Vv these by a constant coupling and constant 
gluon mass of expected size fiy ~ Aqc d ■ In perturbative QCD the effective coupling runs as the inverse of the the 
logarithm. We discuss the effect of this running in the section devoted to the results. Furthermore, form the string 
confinement models it is well known that the string tension changes with the scale as well. At least, it is considerably 
large for the charmonium then for bottomonium. Here we found advantageous to consider such running even for 
different charmonia and we simply separate IS state from the others in the following way 



(5) 



where numerically 6 = 0.9. 

We transform sixteen component BSE into the coupled set of four integral equations for four components A, B, C, D. 
After the projection by jTrj^ we get in Minkowski space : 



A{q,P) 



d*k 



(2^)4 

Sv{-)Sv{+) 



Ka{Vs~Wv) 



-A{k,P) 



p2 

4 



Ssi-)Ssi+)A{k, P) + Sv{-)Ss{+) k.PB{k, P) 



2D{k,P) {k^P'^ + {k.P f 
-C{k,P) 



p2 

2 



(6) 



where we have introduced shorthand notation S{±) — S{k ± P/2). Projecting BSE by Tr^ one gets 



q^B{q,P)+q.PC{q,P) 



_cPk_ 



iq.PKBl + q.kKB2) {Vs + 2Vv) 



Kbi = S^{-)S^{+) 

Kb2 = -SA-)s,(+) 



k P 

-fsik^p) 



k^ + ^]C(k,P 



p^ 



2k.PC{k, P)+[k^ + — S(fc, P 



+ Ss{-)Ss{+)C{k, P) - AeS^{-)Ss{+)D{k, P) 

h Ss{-)Ss{+)B{k, P) + 4k.PS4^)Ss{+)D{k, P) (7) 



Projecting BSE by Tr /P754 one gets in Minkowski space 

74 1 



q.PB(q,P) + C{q,P)P'=t J -^KciVi 



-2Vv) 



Kc 



-S,{-)S,{+)B{k, P)k.P{e - P'^IA) + SsHS,{+)B{k, P)k.P 

p2 



s4-)S4+)Cik,p) 



2{k.py - k^p^ - 



+ Ss{-)Ss{+)C{k,P)P' - S4-)Ss{+)A{k,P)P' 



(8) 



where we have drop out the trivial term proportional to D due to the identity S'y(— — Sv{+)Svi—) — 0. 
Finally the "equation for D" reads 

4 [{q.Pf - g2p2] ^(^^p) ^ .0_K^^P.qP.kVs - P^{k.q)Vs) 
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Kd = Sv{-)Sv{+) 



A{k,P) 



{k^ - P'^/4)D{k,P) 



- Sv{-)Ss{+)C{k, P) + Ss{-)Ss{+)D{k, P) 



(9) 



Using a single homogeneous BSE one automatically disregard resonance character of resonances, i.e. the width 
and presence of the other resonances (couple channel effects) are not taken into account. In this case a quarkonium 
spectroscopy is given by a discrete set of timehke four-momenta — for which BSE has canonicahy normahzed 
solution. To find the solution numerically, it is convenient to perform Wick rotation in the relative momenta, while 
keeping P^ real and positive. In the rest frame of quarkonium P — {M, 0, 0, 0) and after the continuation 



we get for the "Euchdean" Eq. 

Aiq,P) = I ^UAG,iVs-Wv) 
Ua = A{k,P) (k 



+ — +ml\+ 2D{k, P) (-klP"^ + {kiMfj + mck4M(3{k, P) - m^— C(fc, P) , (10) 

where A[q^ P), Ua, ■■■ are the functions of Euclidean momentum qe = {"i-qo, q) = (^4, q) and Minkowski momentum P. 
From now we restrict to the use of free quark propagator 



and introduce shorthand notation for the product of the propagators 
Thus we can simply write 

Ss{-)Ss{+) = mlG2 S,{~)Ss{+) = m,G2 , 



(11) 
(12) 

(13) 



which identities we will use from the next. 

It is easy to find that G2 is the product of two complex scalar propagators , which for equal mass case, becomes 
purely real 



G2 = 



(4 - — + mlf + kjP^ + A„/ 



(14) 



where we have introduced infrared regulator Xinf << rric-, M for a later numerical purpose. Furthermore, we have 
introduced new real function /3 as 



(15) 



which refiects the fact that the function B is odd in the variable k.P — ik^P (valid for positive C parity eigenstate).The 
functions A, /3, C and D are then manifestly real. For purpose of completeness we write the Euclidean Eqs. (|6I7I8I9P 
here 



qlP{q,P)+qiMG{q,P) 



Ub 



{qiMUBl + qE-kEUB2) {Vs + 2Vv)G2 



kAM 



I3{k, P) + y-k% + — + Kj G{k, P) + ^mck%D{k, P) 
Ub2 = {k%- P^/i + ml) P{k,P) + 2kiMG{k,P)-4:kiMmcD{k,P) 



(16) 



M(g,P) + C(g,P)M] 



-Uc{Vs + 2Vv)G2 



Uc = 



p2 

ki{kl + — + ml)p{k, P) + M 



2ki 



p2 

4 



k% + ^+m'^ 



C{k,P) - MmcA{k,P) (17) 



[~iqiM)^ + qlM^]D{q,P) - J {-M^q^kiVs + P\kE.qE)Vs) G2 



Ud = -A(fc,P)-m,C(fc,P) + (-fc|-P74 + m2)i?(fc,P) 



(18) 
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In order to study excited states we keep full k.P dependence into account and do not perform any Sdimcnsional 
reduction which can scrutinize a correct ofF-shell behavior of the vertex function. For given bound state, characterized 
by the mass M the four functions A..D depend only on the scalars q.P and respectively, the first product mix 

Minkowski and Euclidean metrics and become complex {q.P > iq4M in CMS). Due to this fact we we prefer 

to leave 94 integral variables separated from the spacelike part of the integration momentum. When reducing 4-dim 
integral equation we therefore use the following angle integrals: 



J-i [kE - qE) 



{kE - qsY + 

If = f ^^ ,,^'"-^"r\i. ' (19) 

^ y_i [{kE - qEf + IJ^f 

k a 

where z is cosine of angle between space product of momenta z = j^q^) k.q = k^q^ + k.q. We need N = 1,2 ior our 
purpose, for which cases the integrals read 

/W - 1 iJ+ ■ 

- 2\ks\\qs\ ^-' 



l2 , „2 , ,,2 



where 

^± = fc^ + g2 + ^2 _ 2^^^^ ± 2 |k| |q| . (20) 
Substituting Vs,v and integrating over the variable z we get 

dk^ [qiMUBi (k4'' + ^9^Iv) + Ub2 (k4'' + ^9^lP)\ G2 (22) 

(^^^t^('^4"+V4")G2 (23) 



/OO /-OO Ji_"|.2 

(iA;4 / ^C/fl(-MVfc4K4''+7''«4'')t?2. (24) 
-00 Jo v-^^j 



III. NUMERICAL METHOD AND RESULTS 



Contrary to Quantum-mechanical approach, in Quantum Field theoretical framework two different excitations are 
not described by an orthogonal BS functions. Also the number of nodes in various component of BS wave function 
is not driven by any obvious rule. Thus, for instance the dominant component A{p,P) remains node-less for the 
all observed excited states. As different components of BS amplitude vary differently on P.q variable, we do not 
explore more or less conventional expansion into the orthogonal polynomials, which loses its efficiency when, as one 
expects, relatively large number of polynomials could be necessary to distinguish correctly between two excited states. 
Therefore instead of this, we rather solve the full two dimensional integral equation by the method of simple iterations. 

For the purpose of numerical solution we discretized and step by step we scan P^ region of the total momenta. 
Within the step of few MeV we are looking for the solution of the BSE with given P^. Performing several hundred 
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iterations for each given value of we identify the solutions as those for which difference between two consecutive 
iterations vanishes. 

The BSE for bound states is a homogeneous integral equation and it satisfies canonical normalization condition: 

rf^Tr^ J 0^T-,{k,Q)Si+)r,{k,Q)Si-) + ... (25) 

where = and the trace is taken over the Dirac matrices. The conjug ated vertex Tp{k, Q) = C-i^T{-p, Q)C-^ 
with the charge conjugation operator C and three dots represent terms with derivatives of the kernel with respect to 
P . The normalization condition (1251) must be necessarily considered when bound state transitions are considered. 
On the other hand the efficiency of numerical procedure is enforced by the use of the auxiliary normalization which 
is called at each step during the iteration process. For this purpose we implement the auxiliary function \{P) and 
solve numerically the equations (j2ip with A(P) implemented in. For further details of our numerics see Appendix A. 
The physical normalization condition can be applied afterwords. 

Independently on the value of the couplings, we are facing to a rather intriguing problem of BSE equation, which 
we call a doubling of the spectrum of excited states. With the exception of the ground state (and possibly 2s state) 
there are two times more excitations than in the nonrelativistic limit of given model (for nonrelativistic limit see for 
instance [H, [25l - [27| ). We have found that mentioned doubling is numerically stable and we argue that this a rigid 
property of BSE numerical solution for heavy quarkonia within the use of free quark propagators approximation. 
After presentation of the numerical result we come back to this question in the next section , where a simplified model 
with dressed quark propagator will be solved explicitly. 

Nowadays, up to the first two, the excited pseudoscalars have not been accessible in the experiments. Remind 
here, the highest mass candidate suggested is r]{2,s) observed in the Belle experiment in 2007 and hopefully will 
be confirmed in future experiments in recently built FAIR facilities. Nevertheless we do not expect doubling of 
nonrelativistic states in future observed pseudoscalar charmonia. It is well known that the masses of ^ excited 
states should be approximately degenerated with S — states vectors charmonia. The five or six s-wave dominated tp 



charmonium vectors are recently known with masses spread in the energy range 3 — 4.5(5)06^ 5lJ. No doubling is 
observed there and thus we do not expect doubling in the other channels as well. 

In the following we concern the doubling phenomena and explain the identification of the energy pairs. The 
appropriate numerical procedure is described in the Appendix A with a typical numerical search shown in Fig. I VII 
Whenever the "eigenvalue" A cross the unit we get a bound state. The doublets can be characterized by an almost 
identical vertex functions as they are physically distinguishable only by different M . Numerically we found them as 
neighbors with the opposite derivative dX/dM. 

Remind here, the BSE vertex function depends on two variables -.P.q and q^. Alternatively and conventionally, 
vertex function can be visualize as a two dimensional manifold above q^ and q plane. The example of BSE solution 
for ri{lS) is shown at fig. 1 for a fixed p4 = 0.25GeV. The doublets have been actually searched by by comparison of 
the vertex functions. The examples of few neighbor excitations are shown at Figures 2,3,4. The order of the energy 
levels are used to label the solutions (from low to high). 

For a rough estimate we adjusted charm quark pole mass to be ttIc — 1.5GeV^ and we were searching for the 
charmonium spectra for different parameters k{P), /i. The strength of the effective vectorial interaction as = /{Air) 
has been adjusted in order to adjust the intercept to the experimentally known value 1^(28) — rjdls) = 660GeV. 
According to the nonrelativistic calculations we were varying the parameters k between 2 — AGeV^ and /i — 0.1 — 
0.5GeV. At the end the the numerical data have been rescaled in order to get correct value of the ground state 
exactly. The resulting search is presented in the Table 1. , where the all pairs belonging to the single nonrelativistic 
counterpartner are matched. Clearly, all the excited states and even the ground state are heavier then the sum of 
constituent quark antiquark pair in this model. 

Another numerical observation is that the vectorial interaction must be largely suppressed in our model which gives 
typical value estimate Us — 0.07 — 0.1. This is several times smaller then the expected value of the strong coupling 
in nonrelativistic models, where typically one gets as{mc) = 0.5. It is difficult to trace the source of such large 
suppression, it can be due to the sensitivity of BSE to the running of various coupling i.e. k{M) and as{M). 

Here, just in order to exhibit the BSE sensitivity to the running of as we solved the BSE with the running strong 
coupling. For this purpose we use the following simplified separable form 

as{q^k-p):^as{k,p) = — ^''^2 ,„2 n ■ (26) 

In e + -ft-^ 

Within the coupling (pS)) implemented into the BSE with the parameters as(0) = 0.4, Aqcd = 350Mel^ we got 
almost identical data as in the constant coupling case. Clearly, as we reproduce as(Af^^) again, it does not solve 
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Bethe-Salpeter components of r| (IS) 
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FIG. 2: Dependence of components B, (7, D on k of two excited states shown for two slices with given ki — 2.6Me V (I) and 
fc4 = 0.25GeV (II). 



the problem, but it shows up a large sensitivity of BSE to the analytical structure of the kernel. More complete 
understanding of this problem remains for the future numerical study. 

For a ground state mesons it is well known the function A largely dominates. We have found this is the case for 
excited mesons as well. This fact is apparent from the Fig. IVII ?. where the numerical results are compared. 
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Dependence of A on relative "energy" 
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FIG. 3: Dependence of dominant component A of three neighbor excited states shown for fixed threemomenta k = Q.2MeV . 
Vertex functions become wider for higher mass of r/c meson. Quarks are more off-shell in high excited mesons. 



IV. BSE MINKOWSKI SPACE SOLUTION WITH CONFINING PROPAGATORS 

The doubling phenomena is not supported by the experimental data, therefore it should be an artefact of the 
approximations made. Actually, since the BSE model is based on reliable Poincare generalization of the successful 
quantum mechanical models and since the phenomena is qualitatively independent on the interaction details, it 
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TABLE I: Up to date the best fit of conventional BSE solutions for -qdnS). We use conventional quantum mechanical assignment 
nS in order to label states that we expect in nonrelativistic or "instantaneous" approximations. The first column represent 
actual numerical solution in units where rUc = 1.5Gel/, k = 2.%AQGeV^ , in the second column the experimental value of r]{lS) 
has been used to scale other levels. The doubling appears for the states n > 2, and the energy doublets are identified by 
comparison of vertex functions (e.g. by number of nodes in B,C,D functions). After rescaling we produce experimentally 
known r?(25'). For other states -to make levels meaningfully comparable with quantum mechanical labeling- the masses of 
energy levels are averaged for given energy doublets, as = 0.07407 . *Belle observed X(3940) in e^e~ J/ip + X, for the 
interpretation see p^ . 
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FIG. 4: Dependence of B and C components of three excited states shown for fixed threemomenta k = 0.2MeV. Levels labeled 
by B(C)5 (fifth energy excited state) and B(C)6 (sixth excited state) belong to the neigbouhred levels, they have approximately 
identical vertex functions. They have same nodes (odd function B has three nodes, while C functions have two nodes and zero 
minimum at beginning, next two higher levels have only one zero, only seventh is is shown). 



strongly suggests that the doubling effect in calculated quarkonium spectrum is the inherent property of the BSE 
within the use of free propagator approximation. If it so, then the doubling phenomenon should disappear when a 
dressed propagators of confined quarks are taken into account. In this Section we support the idea by an explicit 
numerical solution of very simple "confining" quarkonium model. 

We assume, in opposite to unconfining case, that the propagator of confined object has not a real pole corresponding 
to the free particle solution. Instead of, the singularities of such propagator move aside from the real axis of the 
complex momentum plane due to the strong self-interaction (for example of the electron propagator in QED^ see 
[4l[). The so-called Stingl propagator [13,1131, which exhibits a pair of complex conjugated poles, was suggested for 
the description of a short time living confined gluonic excitations. Such gluon propagator violates reflection positivity 
as it does not satisfy Kallen-Lehmann representation, hence the gluons do not realize in physical spectrum. Actually, 
more recently (4^ . l45j it was discussed how such propagator can quantitatively describe lattice predictions for the 
infrared Landau gauge gluon propagator [46l |47|. Remind also here, that the structure of the quark propagator pro 
has been estimated in Euclidean approximation (48j . where the quark propagator calculated by Schwinger-Dyson 
equations has been fitted by the sum of several complex conjugated poles. 

Here we extend Stingl's idea to the case of quarks as well and approximate quark propagator by the following 
formula 



(fc2 - m2)2 + 54 



(27) 



which is written in Minkowski space and where rric, S are real parameters and / is a complex constant. For quark with 
momentum k ~ (5 the parameter 6 can be interpreted as minimal wavelength [52| . In general / is complex function of 
the momenta, which in accordance with perturbation theory must become a real valued function for a large spacelike 
k. Here, however we simplify and take a constant phase approximation 



/ = exp^ 2 



(28) 



in order to approximate its infrared behavior only. Dealing with heavy quarks, we also assume S is smaller then rric, 
rric — l.hGeV in accordance with other QCD estimates. 

The BSE is usually solved by performing the so called Wick rotation at the first step and the vertex functions are 
calculated for Euclidean relative momenta. This is well defined treatment when the kernels and the functions we are 
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FIG. 5: Dependence of dominant component A of two excited states shown for two fixed values k4 = 2.6AleV and k4 = 0.25GeV . 
The function A has no nodes for Euclidean relative momentum k. 



looking for are analytical inside the integration contour. It is certainly not the case of many models considered in 
the literature and Minkowski solution is sometimes defined as an analytical continuation of the Euclidean one, i.e. 
Euclidean space can be regarded as a starting point albeit the analytically continued solution can differ from the 
originally defined Minkowski one. The propagator with complex conjugated poles is just an example. 

We avoid a problematic Euclidean to Minkowski space continuation of numerical BSE solution and solve the 
BSE directly in Minkowski space. It is worthwhile meaning, that this is well defined numerical problem and no 
troubles appear, which would eventually prevent the numerical solution (note here, that when trying to solve a sort 
of Schwinger-Dyson equations in Minkowski space, then apart an usual problems related with various pole singularity 
structure, one is also faced to hard technical difficulties related with improper principal value momentum integral, 
this trouble does not arise here, because of asymptotic of the bound states vertex function) . 

To make comparison meaningful we keep the Lorentz structure of the kernel in the form as in the previous Euclidean 
study. However to make the model more easily tractable we change the double pole of the scalar kernel to the complex 
conjugated one as well, while we leave the vectorial in its original single pole form. More explicitly we take now: 



K = , (29) 

The detailed derivation will be published elsewhere for its own interest, here we stress the main difference and 
results. We stay in Minkowski space, the square of momenta k reads fc^ = /cq — k^j and the BSE turns to be two- 
dimensional integral equation for two real variables: two scalars k^,k.P for given discrete = or alternatively 
for the off-shell relative energy /cq and spacelike relative momentum k. As in the previous case we assume that the 
use of the function A is enough for a reasonable estimate of the spectrum and we solve the BSE for this component 
in very similar fashion as in the previous case. 

Numerically we take now /Xs = /i„ = /i = 0.256*6^ and 2/i for our first numerical estimate. Having double 

pole changed to the the complex conjugated poles the strength is effectively suppress. We found that taking As too 
small can lead to numerical instability, while taking As too large decreases the coupling strength and population 
of energy levels as well. To get a reasonable solution of BSE we increase the scalar coupling strength, numerically 
we take k = 5GeV^. Within rric — l.5GeV the spectrum we have found reads : M(ls) = 1.97GeV; M{2s) — 
3M3GeV;M{3s) = i.WSGeV; M{As) = 5.9GeV; M{5s) = 6MGeV. The doubling in the numerical solution of BSE 
disappeared. The intercept is presented, although it is overestimated. This quantitative disagreement calls for an 
improvement of the approximation made and further possible refinement of the presented model. 
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V. CONCLUSION 



The homogeneous BSE for excited states of pseudoscalar charmonia have been solved. When we use free -unconfined- 
approximation for the quark propagators, we get numerical solution for arbitrarily high excitation, however when 
comparing with nonrelativistic solutions the excited states are overpopulated by the factor 2. The excited states come 
in pairs with close masses and almost identical wave functions. This doubling phenomena is very likely an artefact of 
the approximations. In the model with confining form of propagators the doubling of the energy levels disappeared. 
However as the model suffers by further simplifications and approximations it does not reproduce the experimental 
data in recent stage. Nevertheless, it is reliable enough in oder to illustrate its main effect- disappearance of the 
energy levels doubling. 

Since the high number of excited states has been found by the numerical solution of BSE method and the author is 
not aware about similar numerical findings in the literature, the numerics has been described in details in the paper. 
For those who are interested the code is available by an email request. 



VI. APPENDIX A 

In this Appendix we write down the BSE in the form which has been actually used in our numerical solution. We 
found it was suited to rewrite Eqs. (|2m24l) in the following equivalent form: 

/oo />oo Juu2 . , 

[-((74M)2 + g|M2]Vei-oo Jo (2^)' ^ ^ 

l3{q,P) = A(P)— ^%^^l-!-^^^(lhs. 0fEq.)g2 + (rhs. of Eq. i22))q,M 

C(q,P) = \{P) "(g4^^) +ggM ^^^^^ of Eq.)p2 + (rhs. of Eq. {2i))qiM , (30) 

[-{qiMf + q^AP] + e 

where e is a small numerical regulator and where the last two equations are equivalent to the Eq. (22) and Eq. (23) 
in the limit e — > 0. The integrals have been discretized and the system of equations solved by iterations as usually: 
The left hand side of Eqs. represents i — th iteration step result, when i — 1 — th iteration has been used to evaluate 
Ihs. of Eqs. system. Starting with reasonable zeroth approximation then 300-400 iterations have been used for each 
fixed value of P^, which was enough to get obviously convergent solution. 

Main goal is the implementation of the function X{P) which when properly chosen can radically increase the 
efficiency of the numerical search. Among several working possibilities that we have checked we mention two: firstly 
the choice A^^(P) = ^(0,0) appears to work well in many cases. As the second, we have taken 



which has been actually used for the evaluation of the results in the presented paper. The solution of the BSE are 
identified whenever A = 1. Non-trivially, when A = 1, the difference between consecutive iterations vanishes at the 
same time. 

The example of the solution of BSE is shown in Fig lVIl . The function sigma is the weighted integral difference 
between two consecutive iterations evaluated through the following formula: 

To get a reasonable numerical error we found that relatively large number of integration points is required. Therefore 
to speed up numerics we identify roughly the positions of the bound states within small -40 * 80- number of mesh- 
points. Consequently the results has been improved within -88 * 176- points of discretized integral variables k, fc4 
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FIG. 6: Numerical search for = l.5GeV, k = 2.5GeV'^,as = 0.08, /x„ = /^s = 350Mel/ (solid line). Dashed (dot-dashed) 
line represents A, (a) in the approximation when only A vertex function has been considered, in this case rric — 1.5GeV, k = 
2.5GeV^,as = 0.0815, ^„ ^ fJ-s ^ 350MeV. 
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